PLOS Computational Biology

Detailed agent-based simulations are increasingly used to support policy decisions, but their computational cost and complex uncertainty structure make systematic scenario analysis challenging. We present a data-driven, uncertainty-aware decision support (DDUADS) workflow for using stochastic simulation models as decision-support tools under limited computational budgets. The approach combines several established techniques-sensitivity screening, Bayesian calibration using simulation-based inference, and multi-surrogate model integration for translational efficiency-into a coherent pipeline that enables uncertainty-aware policy analysis. Rather than producing a single baseline, the calibration stage yields a posterior distribution over plausible model parameterizations, allowing flexible, uncertainty-aware forward projections. We demonstrate the DDUADS workflow on the INFORM-HIV agent-based model of HIV transmission in Chicago to evaluate potential disruptions in antiretroviral therapy (ART) and pre-exposure prophylaxis (PrEP) use. While the specific application is HIV modeling, the challenges and techniques described here arise in other simulation studies and can be applied to decision support in other domains.

Patient-specific four-chamber electromechanical models provide a physics-constrained framework for investigating whole-heart cardiac physiology and disease mechanisms. Identifying which model parameters impact whole-heart function is important for understanding cellular-, tissue-, and organ-scale determinants of cardiac performance and for calibrating patient-specific models. However, previous global sensitivity analyses of cardiac electromechanical models have typically been performed on a single heart, and systematic evaluation of how parameter influence compares across anatomically different subjects remains limited. We created four-chamber electromechanical models using cardiac MRI from five healthy subjects (n = 5). The models simulated atrial and ventricular cellular electrophysiology, calcium dynamics, and active contraction, with heterogeneous fibre orientation, transversely isotropic tissue mechanics, pericardial constraint, and a closed-loop cardiovascular system providing physiological boundary conditions. In total, 46 parameters described the integrated model. Using Gaussian process emulators, we performed multi-scale global sensitivity analysis to evaluate the relative contribution of model parameters to left and right atrial and ventricular function. Across all anatomies, the most influential parameters were systemic and pulmonary resistances, ventricular end-diastolic pressures, and the venous reference pressure, highlighting the dominant role of haemodynamic loading conditions in governing pressure- and volume-based outputs. A chamber-level analysis of atrioventricular coupling revealed a phase-dependent pattern. Atrial pressures were predominantly governed by global haemodynamic parameters (> 90% of total sensitivity), atrial filling volumes showed substantial ventricular influence ({approx}40-55% across anatomies), and atrial end-systolic volumes were primarily determined by intrinsic atrial parameters ({approx}60-65%). These patterns were consistent across subjects despite differences in anatomy. We show that, in healthy male subjects, inter-individual anatomical variation does not substantially change the ranking of dominant parameters. This work provides a repeatable modelling and sensitivity analysis framework and establishes a benchmark reference for whole-heart electromechanical modelling in healthy hearts. Author summaryComputational models of the heart can simulate cardiac physiology in unprecedented detail, but these models contain many parameters whose influence on predicted function is not fully understood. We built patient-specific four-chamber heart models from MRI scans of five healthy subjects and used statistical methods to systematically test how 46 model parameters affect simulated cardiac performance. Across all five subjects, we found that the haemodynamic loading parameters, including systemic and pulmonary vascular resistance, ventricular filling pressures, and the venous reference pressure, consistently had the greatest influence on the model outputs, regardless of differences in individual heart anatomy. This finding suggests that in healthy resting conditions, the boundary conditions of the cardiovascular system, rather than individual differences in heart geometry or electrical properties, are the primary drivers of whole-heart function. We also found a structured coupling pattern between the upper and lower heart chambers, where global haemodynamic parameters dominate atrial pressure regulation, ventricular mechanics shape atrial filling, and intrinsic atrial properties control atrial emptying. This work provides a benchmark dataset of five anatomically detailed heart models and a sensitivity analysis framework to guide calibration of future cardiac digital twin models.

Neuronal activity alters extracellular ion concentrations and electric potentials. Ephaptic effects refer to the feedback influence that these extracellular changes can have on neuronal activity. While electric ephaptic effects occur on a fast timescale due to extracellular potential perturbations, ionic ephaptic effects are driven by slower, accumulative changes in ion concentrations. Among the previous computational studies of ephaptic effects, the vast majority have focused exclusively on electric effects, while ionic ephaptic effects have largely been neglected. In this work, we present an electrodiffusive computational framework consisting of two-compartment neurons that interact via a shared extracellular space. By accounting for both electric potentials and ion-concentration dynamics in a self-consistent manner, our framework enables us to explore the relative roles of electric and ionic ephaptic effects. Through numerical experiments, we demonstrate that ionic and electric ephaptic interactions play very different roles. While ionic ephaptic interactions increase population firing rates, electric ephaptic interactions primarily drive subtle shifts in spike timing. Furthermore, we show that these spike shifts cause the phase difference (the distance in spike times between a small collection of neurons) to converge to a stable, unique phase difference, which we coin the ephaptic intrinsic phase preference. Author summaryNeurons predominantly communicate through synapses: specialized contact points where a brief electrical signal, known as a spike or action potential, in one neuron influences another. Neurons generate these spikes by exchanging ions with the surrounding extracellular space. This way, spiking neurons alter extracellular ion concentrations and electric potentials. Since neurons are sensitive to such changes in their environment, they can also influence one another indirectly through the shared extracellular medium. This form of non-synaptic interaction is known as ephaptic coupling. Most computational models of neuronal activity neglect ephaptic interactions, and those that include them typically consider only electric effects while ignoring ionic contributions. As a result, the relative roles of electric and ionic ephaptic effects remain poorly understood. Here, we introduce a computational framework that accounts for both mechanisms in a self-consistent way. Our results show a functional distinction: ionic ephaptic effects act slowly, regulating population firing rates, whereas electric ephaptic effects act on millisecond timescales and subtly shift spike timing. These shifts cause spike-time differences between neurons to converge to a stable value, a phenomenon we call ephaptic intrinsic phase preference.

Genetic surveillance is increasingly used to track malaria transmission, yet genomic metrics can respond nonlinearly to changes in transmission intensity and depend on the diversity already present in the parasite population. Here, we present a stochastic agent-based model of hu-man-mosquito transmission that integrates SEIS-like epidemiological dynamics with within-host Plasmodium falciparum haplotype dynamics. By varying the maximum mosquito biting rate and the initial parasite diversity, we examine how transmission intensity and standing diversity jointly shape mixed infections, recombination, and long-term population structure across a continuous transmission gradient. Our study revealed a sequential pattern in which increasing biting intensity first increases infection prevalence and multiplicity of infection, then expands opportunities for outcrossing, and only thereafter increases effective recombination and recombinant haplotype generation. These responses are strongest in low- to intermediate transmission and tend to plateau at higher transmission levels. Initial population diversity constrains the amount of diversity that can be maintained and the magnitude of recombination output, while temporal trajectories show that haplotype evenness can pass through transient non-equilibrium phases before stabilizing. Together, these results show that the structure of the parasite population is shaped not by trans-mission intensity alone but by its interaction with standing genetic diversity. Furthermore, this study works to clarify when and how genomic metrics reliably reflect transmission conditions across heterogeneous malaria settings.

Evidence accumulation is a fundamental aspect of human decision-making. However, how the precise temporal structure of evidence shapes the accumulation process has not been systematically studied. As a result, current understanding of evidence accumulation remains largely limited to its time-averaged behavior. We tested human subjects in a visual estimation task in which they inferred the angular position of an unknown source from a noisy stimulus sequence. Introducing systematic temporal perturbations, i.e., breaks of different durations and at different positions in the otherwise regular evidence sequence, revealed that subjects actively compensated for the memory loss endured during the break by dynamically enhancing evidence integration and memory maintenance immediately after the break. We derived a new time-continuous Bayesian updating model that is dynamically constrained by optimal performance-effort trade-offs. With two free parameters determining the overall resource-efficiencies of encoding and memory maintenance, the model accurately predicts the rich dependencies of subjects accumulation behavior on the evidence schedule, including subjects individual tendencies to emphasize either early (primacy) or late (recency) samples in the evidence sequence. Our results suggest that evidence accumulation is a non-stationary, dynamically controlled process that optimally balances the information gained from incoming evidence against the cognitive effort required to acquire and maintain it. The proposed model is general and should apply broadly across many task domains.

Positive selection on somatic mutations is the driving force for cancer progression. Growing evidence shows that the emergence of a driver mutation in a tumor sample depends on individual-specific factors, for example environmental exposures or the individuals germline genetic background. We term these individual-level factors as the "contexts" of a tumor. Our hypothesis is that mutations in a driver gene can bring different growth advantages in different contexts, resulting in "differential selection" on these genes in varying contexts. Identifying which contexts modulate selection strength provides critical insights into the selection forces driving tumorigenesis. However, due to the sparsity of somatic mutations and heterogeneous background mutational process across positions and individuals, identification of differential selection has limited power with current statistical tools and is prone to false positives. To address this, we developed a powerful statistical method, DiffDriver, that identifies associations between "contexts" and selection strength on a driver gene across individuals. DiffDriver accounts for variations of mutation rates across bases and individuals, while taking advantage of functional information of sequences to improve the power. Through simulations, we show DiffDriver reduces false positives and boosts power compared to current methods. Our results highlight that multiple individual-level factors create significant heterogeneity in the strength of selection acting on driver genes and 33% of driver genes showed differential selection in at least one of the contexts studied, including tumor clinical traits and tumor immune microenvironment subtypes. These results provided new insights into the context-dependent forces driving cancer evolution.

Neural population models are widely used to interpret electroencephalography (EEG), yet their relationships remain far less systematically understood than those among single-neuron models. More fundamentally, it remains unclear whether EEG can support a uniquely plausible population-level mechanism, or whether multiple structurally distinct models can explain the data equally well. To address this question, we combine comparative analysis of canonical model families with grammar-based generation of new candidate architectures. We assembled 17 canonical neural mass and phenomenological models and embedded them in a shared structural space. From their common processes, we defined a probabilistic grammar over interpretable dynamical components and developed ENEEGMA (Exploring Neural EEG Model Architectures), a Julia-based framework for grammar-based model generation, simulation, and parameter optimization, to generate additional candidate models. We then assessed both canonical and generated models by fitting them to EEG independent-component spectra from four datasets for each condition: resting state and steady-state visual evoked potentials. Canonical models formed six structural clusters. Across conditions, compact low-dimensional polynomial oscillators performed best overall, with Montbrio-Pazo-Roxin, FitzHugh-Nagumo, and Stuart-Landau models offering the best balance of fit quality, stability, and simplicity. Grammar-based exploration further showed that the space of viable EEG node models extends beyond canonical formulations: even a restricted search over 1,000 generated models produced compact alternatives competitive with nearly all canonical families and achieving the strongest cluster-level SSVEP fits. Together, these findings suggest that EEG spectra constrain plausible neural population mechanisms without uniquely determining them. Beyond this, grammar-based model exploration provides a principled, data-driven framework for EEG-constrained model discovery. Author summaryElectroencephalography (EEG) lets us measure brain activity non-invasively, but the signals are indirect, so we rely on mathematical models to explain how neural populations generate them. Many such models exist, yet it is unclear whether standard models cover the full range of plausible explanations for EEG data, or whether several very different models can explain the same signal equally well. In this study, we compared a broad set of established neural population models and then used a grammar-based equation discovery framework to automatically generate new candidate models from interpretable building blocks. We found that simple low-dimensional oscillator models often matched EEG spectra better than more complex canonical models. We also found that newly generated models could perform nearly as well as, and sometimes better than, established ones, especially for stimulus-driven responses. These results suggest that EEG spectra alone may not be enough to identify a unique underlying neural mechanism. More broadly, our work shows how automated, biologically informed model generation can help to compare, understand, expand, and test the space of candidate neural population models.

The connectivity of brain networks is extremely sparse due to metabolic, physical and spatial constraints. Although wiring sparsity can confer computational advantages for biological and artificial neural networks, sparse networks require fine parameter tuning and exhibit strong sensitivity to perturbations. How brains achieve their efficiency and robustness is unclear. Here we addressed this by analysing the dynamical properties of Echo State Networks with wiring based on the Drosophila melanogaster fruit fly connectome, compared with sparsity-matched random-wiring networks. We evaluated these networks on a set of eight cognitive tasks, and found that connectome-based neural networks (CoNNs) typically showed narrowly distributed task engagement across their neurons. The importance of a neuron for task performance correlated with its node degree, local clustering, and selfrecurrency, and these correlations were stronger in CoNNs than in random networks. CoNNs were more robust to neuronal loss, retaining their task performance and beneficial dynamical properties such as criticality and spectral radius better than random networks. Similarly, CoNNs were more robust to hyperparameter variations in both input and recurrent weight scaling. Using theoretical arguments and numerical simulations, we show that excess CoNN node self-recurrency is sufficient to explain this enhanced robustness. Overall, these results identify non-random features of connectome wiring that allow brains to reconcile extreme sparsity with reliable computation. SignificanceBrain networks support robust computation even though they operate under extreme wiring sparsity due to metabolic and spatial constraints. While sparse networks typically require fine-tuning and are sensitive to perturbations, we show that biological connectomes support specialised, efficient task engagement and remain robust to neuron loss and parameter variation. We identify excess neuronal selfrecurrency as a key structural feature underlying this stability. These results reveal how non-random connectivity stabilises computation in extremely sparse networks, providing principles for understanding brain function and designing robust, efficient artificial neural systems.

Developing approaches to link structure and function is an ongoing challenge in computational and structural biology. Using a systems-level framework, we present here an analysis pipeline in a Python package, mdsa-tools, that constructs network representations of structures in a time series of trajectory frames from molecular dynamics (MD) simulations. Here, we demonstrate its use on a ribosomal subsystem. The subsystem is centered on the CAR interaction surface, a "brake pad" adjacent to the aminoacyl (A-site) decoding center that tunes protein translation rates. We leverage unsupervised learning algorithms to explore the conformational landscape of behaviors visited by two versions of the subsystem (brake-on and brake-off) that differ at the codon 3 adjacent to the A-site codon. Our network representations of MD frames embody H-bond interactions between all pairwise combinations of residues in the system. By utilizing per-frame vector representations of network edges, we can apply standard clustering and dimensionality reduction methods to explore behavioral differences between the brake-on and brake-off versions of the system. K-means clustering of frame vectors revealed a striking separation of the two system versions, consistent with principal components analysis (PCA) embeddings and Uniform Manifold Approximation and Projection (UMAP) embeddings. Dissection of K-means centroids and PCA loadings highlighted H-bond interactions between residue pairs in the ribosomes peptidyl site (P site), suggesting potential allosteric signaling across the subsystem. Author summaryWith the impressive development of computational algorithms to successfully simulate the dynamics of biological molecules over time, the exploration and incorporation of systems modes of analysis is a natural next step to begin to understand the molecular dynamics behaviors that emerge from these experiments. Following the approaches of classical molecular genetics, we used a "computational genetics" paradigm where we introduced changes (mutations) in potentially important residues, changing their identities or modifying their chemical properties, and asked how the dynamic system responded to these changes, viewing the simulations as a series of movie frames of the dynamic structure over time. Starting with network representations of each frames structure, where the nodes are residues, and the edges denote H-bond interactions between the residues, we used several unsupervised machine learning algorithms to uncover behavioral changes in the different mutated versions of the system. Applied to our ribosome neighborhood, this revealed unexpected changes in behavior at the ribosome peptidyl site (P site) in response to mutating mRNA residues on the other side of the aminoacyl site (A site) codon, suggesting long-range allosteric interactions across the neighborhood.

Connectome-constrained neural network models aim to link anatomical connectivity with functional computation by training networks whose architectures reflect biological circuits. Because such models are increasingly used to infer neural mechanisms, it is important to assess their robustness to variations in training conditions and model selection criteria. Here we retrain ensembles of connectome-constrained models under nominally identical conditions and compare their correspondence to experimentally measured response properties in the Drosophila motion pathway. While task performance remains similar across models, the identification of biologically plausible circuit solutions is unstable across retraining runs. In particular, model clusters selected by lowest validation task error do not reliably correspond to experimentally observed neural tuning, and small variations in performance metrics can reorder cluster rankings. These results indicate that, in this framework, similar task performance does not reliably identify biologically plausible circuit solutions. Task error alone is therefore insufficient for mechanistic identification, and additional model-selection criteria are needed.

BackgroundFlexibility and robustness of neuronal function are closely linked to degeneracy, the ability of distinct structural or parametric configurations to produce similar functional outcomes. At the cellular level, this often manifests as ion-channel degeneracy, in which multiple combinations of intrinsic conductances yield comparable electrophysiological phenotypes. MethodologyWe used a population-based, data-driven modelling framework to generate large ensembles of biophysically detailed CA1 pyramidal neuron models constrained by somatic electrophysiological features extracted from patch-clamp recordings in acute slices from early-birth rats. 10 reconstructed morphologies were incorporated, and model populations were analyzed using parameter correlation analysis, principal component analysis, and generalization tests to assess robustness, degeneracy, and morphology dependence of intrinsic properties. ConclusionsAcross the model population, similar somatic firing behaviours emerged from widely different combinations of intrinsic parameters, demonstrating robust two-level ion channel degeneracy both within and across morphologies. Each morphology occupied a distinct region of parameter space, indicating morphology-specific compensatory effects, while weak pairwise parameter correlations suggested distributed compensation rather than tight parameter dependencies. Even with a fixed morphology, multiple parameter subspaces supported comparable electrophysiological phenotypes. Generalization across morphologies was structure-dependent and non-reciprocal, with successful parameter similarity occurring preferentially between structurally similar neurons. Interestingly, to accurately simulate spike-frequency adaptation, it was important to retain some kinetic properties of the ion channel models as free parameters during optimization. Together, these findings show that dendrite morphology shapes the valid parameter space, and similar electrophysiology of CA1 pyramidal neurons arises from the interplay between structural variability and ion-channel diversity. This work highlights the importance of population-based modelling for capturing biological variability and provides insights into how neuronal robustness might be maintained despite substantial heterogeneity, and offers a scalable pipeline for generating biophysically realistic CA1 neuron populations for use in network simulations. Author summaryNeurons must reliably process information even though their internal components, such as ion channels and cellular shape, can vary widely from cell to cell. How stable behaviour emerges from such variability is a fundamental question in neuroscience. In this study, we explored this problem using detailed computer models of early-birth rat hippocampal CA1 pyramidal neurons, a cell type that plays a central role in learning and memory. Instead of building a single "average" neuron model, we created large populations of models that all reproduced key experimental recordings but differed in their internal parameters. We found that neurons with different shapes and different combinations of ion channels could nevertheless generate similar electrical activity. This phenomenon, known as ion channel degeneracy, allows neurons to remain functional despite biological variability or perturbations. Our results show that neuronal shape strongly influences which parameter combinations are viable, but that multiple solutions exist even for the same morphology. The population of models we provide offers a resource for future studies of early-birth CA1 pyramidal cell function and dysfunction.

Resource sharing is a fundamental form of social exchange underlying the formation and maintenance of social bonds in humans and other species. While reciprocity has long been proposed as a key mechanism in group interactions, the dynamic processes underlying resource allocation remain poorly understood. In this study, we employed computational modeling to investigate the temporal dynamics of resource sharing in a novel group decision-making task across three experiments. We found that, beyond the well-documented reciprocity, participants exhibited consistent alternating behavior, characterized by the switching between potential recipients. This alternation was not driven by fairness concerns but reflected a strategic balance between maintaining stable partnerships and exploring alternatives. Crucially, a reinforcement learning model incorporating Theory of Mind (ToM) consistently outperformed all alternative models. These findings highlight the critical role of ToM in social decision-making and suggest that mentalizing others intentions may be essential for effective resource sharing and social bond formation.

Coral reef management under climate change is challenging due to data sparsity and high uncertainty, yet it is essential for informing conservation strategies. We present CoralBlox, a mechanistic discrete time coral ecology model with the explicit aim of supporting rapid scenario exploration and decision making. The model represents discretized distributions of five coral functional groups across configurable spatial scales while incorporating key ecological processes, including coral growth, reproduction, thermal adaptation, and responses to disturbances. Validation against observed data demonstrates that CoralBlox effectively captures major trends in coral cover dynamics across the Great Barrier Reef, particularly for bleaching-driven mortality and recovery patterns. While simplifying ecological complexities, the model maintains sufficient ecological realism to evaluate and compare the result of distinct management strategies. CoralBlox enables comprehensive assessment of potential management interventions with high computational efficiency and interoperability. The models flexible architecture makes it extensible to coral ecosystems worldwide, providing valuable exploratory capability for reef management. TeaserCoralBlox is an efficient coral reef ecology model supporting rapid scenario testing and management decision making under climate change. HighlightsO_LIMarine ecosystems are characterized by high uncertainty and data sparsity. C_LIO_LIManagement decisions still need to be made under these uncertain contexts. C_LIO_LICoralBlox offers a conceptually simple yet credible representation of ecological processes. C_LIO_LIComparatively fast runtimes across different spatial scales enable rapid exploration of plausible future states. C_LI

The sodium-potassium ATPase (NKA) consumes 19-28% of cellular ATP and is critical for maintaining ion homeostasis. Understanding its energetic efficiency is essential for comprehending cellular physiology and pathophysiology. We develop bond graph models of the NKA that ensure thermodynamic consistency by enforcing conservation of mass, charge, and energy. A simplified 6-state model captures biophysics comparable to a 15-state model while remaining computationally tractable. Through detailed energetic analysis, we demonstrate that under physiological conditions, approximately 65% of the energy from ATP hydrolysis is stored as chemical energy in ion gradients, 10% as electrical energy in the membrane potential, and 25% is dissipated as heat, yielding an overall efficiency of [~]75%. We investigate how the free energy of ATP hydrolysis ({Delta}GATP), intracellular Na+, and extracellular K+ affect NKA efficiency and activity. A critical threshold exists at {Delta}GATP {approx} - 48 kJ/mol below which chemoelectrical transduction drops dramatically, consistent with NKA inhibition under ischemic conditions. The bond graph framework enables quantitative comparison of different NKA models and provides a systematic approach for analyzing ion pumps. SIGNIFICANCEThe sodium-potassium ATPase is one of the bodys most energy-consuming enzymes, yet its energetic efficiency and mechanisms remain incompletely understood. This study presents the first comprehensive energetic analysis using bond graph modeling, guaranteeing thermodynamic consistency. By demonstrating that simplified 6-state models capture essential energetic behaviors of complex 15-state models, we establish bond graphs as a powerful, tractable tool for energetic analysis, model comparison, model selection and validation. The bond graph approach can be applied to other transporters, offering a powerful tool for systems physiology and drug discovery.

Recent progress in mathematical kinship modelling has allowed one to predict the probable numbers of kin for a typical population member. In the models, kin may be structured by age and sex, both in static or time-variant demographies. Knowing the probable numbers of kin in different stages - such as parity, health status, or geographic location - however, remains an open challenge in Kinship Demography. Knowing how population structure delimits kin to distinct stages is an advance - for instance, the probability of having one sister at home and one sister away has different social implications from the probability of having two sisters. We present a novel analytical framework, grounded in branching process theory, that provides kin-number distributions jointly structured by age and stage. Using recursive compositions of probability generating functions (PGFs), we derive the joint age, stage, and age x stage kin-number distributions. All marginal distributions over either dimension naturally emerge. Simple extensions of the PGF approach additionally yield: the joint distribution of an individuals own stage and their kins stage; the probable numbers of kin deaths, both in total and by generation number; and the probabilities of being kinless and/or orphaned. We demonstrate the framework through novel results in an application using UK parity-specific fertility and mortality data. HighlightsO_LIA new method calculates probability generating functions for the number of kin structured by age and stage C_LIO_LIThe model allows predicting the probable numbers of kin organised by age and stage C_LIO_LIRecursive nesting of probability generating functions in branching processes is used C_LIO_LIAn application is presented highlighting the novel results C_LI

Hidden (semi-) Markov Models (HsMMs) are increasingly being used to segment neurophysiological signals into sequences of latent cognitive processes. The idea: different processes will leave distinct traces in trial-level recordings of (multivariate) neuro-physiological signals. Markov models, equipped with an emission model of these traces and a latent process model describing the progression through the different latent processes involved in a task, can then be used to infer the most likely process for any time-point and trial. However, the currently used HsMMs remain limited in two important ways. First, they cannot account for subject-level heterogeneity in the latent and emission process. Instead, a single group-level model is assumed to explain the entire data. Second, they cannot account for the potentially non-linear effects of experimental covariates on the latent and emission process. To address these problems, we present a modeling framework in which the HsMM parameters of the emission and latent process are replaced with mixed additive models, including smooth functions of experimental covariates and random effects. We derive all necessary quantities for empirical Bayes and fully Bayesian inference for all parameters and provide a Python implementation of all estimation algorithms. To demonstrate the advantages offered by this framework, we apply such a multi-level model to an existing lexical decision dataset. We show that, even in such a simple task, not all subjects rely on the same processes equally and that at least two semi-Markov states, previously believed to reflect distinct processes, might actually relate to the same cognitive process.

Fine-grained temporal structures emerge in smartphone behavioral recordings over multi-day periods. Complex systems research suggests that emergent temporal structures reflect underlying resource constraints of the system. Here we test whether cognitive abilities measured through speeded tasks (spanning fractions of a second) are reflected in emergent smartphone temporal structures spanning days, revealing how cognitive resource limitations shape naturalistic behavior. We analyzed smartphone tap interval patterns accumulated over several days and used decision tree regression models to predict performance in simple and choice reaction time tasks from these patterns. Simple reaction time was poorly predicted (R2 = 0.003), indicating that basic sensorimotor constraints play only a marginal role in shaping real-world behavioral timing. In contrast, choice reaction time was moderately predictable (R2 = 0.4), demonstrating that higher-order cognitive constraints prominently influence naturalistic temporal organization. Notably, while task performance operates at sub-second timescales, predictive temporal patterns in smartphone behavior spanned milliseconds to several seconds and was accumulated over days, revealing the broad, multi-scale influence of cognitive resource constraints on everyday behavior. Both predicted and measured choice reaction times showed age-related decline, but the decline was more pronounced in predicted values, suggesting that age-related cognitive changes may be amplified in naturalistic contexts. These findings demonstrate that emergent temporal structures in smartphone use can reveal how cognitive processes measured using speeded tasks manifest, or fail to manifest, in real-world behavior. These findings demonstrate that complex-systems approaches can bridge laboratory and naturalistic assessments of cognition, revealing which cognitive processes meaningfully constrain real-world behavior.

Contrastive principal component analysis (PCA) methods are effective approaches to dimensionality reduction where variance of a target dataset is maximized while variance of a background dataset is minimized. We previously described how contrastive PCA problems can be written as solutions to generalized eigenvalue problems that maximize particular instantiations of the Rayleigh quotient. Here, we discuss two extensions of contrastive PCA: we use kernel weighting from spatial PCA (k-{rho}PCA) to contrast spatial and non-spatial axes of variation, and separately solve the Rayleigh quotient in the space of basis function coefficients (f-{rho}PCA) to find modes of variation in functional data. Together, these extensions expand the scope of contrastive PCA while unifying disparate fields of spatial and functional methods within a single conceptual and mathematical framework. We showcase the utility of these extensions with several examples drawn from genomics, analyzing gene expression in cancer and immune response to vaccination.

Recent Artificial Intelligence (AI) models have matched or exceeded human experts in several benchmarks of biomedical task performance, but have lagged behind on surgical image-analysis benchmarks. Since surgery requires integrating disparate tasks --- including multimodal data integration, human interaction, and physical effects --- generally-capable AI models could be particularly attractive as a collaborative tool if performance could be improved. On the one hand, the canonical approach of scaling architecture size and training data is attractive, especially since there are millions of hours of surgical video data generated per year. On the other hand, preparing surgical data for AI training requires significantly higher levels of professional expertise, and training on that data requires expensive computational resources. These trade-offs paint an uncertain picture of whether and to-what-extent modern AI could aid surgical practice. In this paper, we explore this question through a case study of surgical tool detection using state-of-the-art AI methods available in 2026. We demonstrate that even with multi-billion parameter models and extensive training, current Vision Language Models fall short in the seemingly simple task of tool detection in neurosurgery. Additionally, we show scaling experiments indicating that increasing model size and training time only leads to diminishing improvements in relevant performance metrics. Thus, our experiments suggest that current models could still face significant obstacles in surgical use cases. Moreover, some obstacles cannot be simply ``scaled away'' with additional compute and persist across diverse model architectures, raising the question of whether data and label availability are the only limiting factors. We discuss the main contributors to these constraints and advance potential solutions.

Local field potentials recorded from the subthalamic nucleus (STN) in Parkinsons disease (PD) exhibit a distinctive multiscale spectral signature: exaggerated beta-band oscillations (13-30 Hz) coupled to high-frequency oscillations (HFOs, 200-400 Hz), with HFO amplitude being phase-locked to the beta cycle. This phase-amplitude coupling (PAC) has been identified as a promising biomarker of the parkinsonian state, yet no biophysical model has explained how it emerges, what determines the HFO frequency, or how HFOs can exist without beta modulation in the medicated STN. Here we show that a heterogeneous population of excitatory Izhikevich neurons with recurrent coupling produces three dynamical regimes: (i) asynchronous tonic firing, (ii) asynchronous bursting, in which neurons burst individually producing broadband HFO power but without coherent population-level PAC, and (iii) synchronous bursting, which gives rise to beta-HFO PAC. The regimes are governed by two biophysically interpretable parameters that capture complementary effects of dopamine depletion: one reflecting changes in intrinsic neuronal excitability, the other reflecting changes in synaptic coupling strength. The transition from asynchronous to synchronous bursting in this model captures the emergence of pathological STN neuronal activity in the parkinsonian state. HFO peak frequency varies continuously across the two-parameter landscape, providing a mechanistic account of the clinically observed shift from slow (200-300 Hz) to fast (300-400 Hz) HFOs between medication states. The character of the synchronization transition depends on baseline excitability, ranging from a sharp co-emergence of bursting and synchrony at low excitability to a decoupled two-stage process at intermediate excitability where burst recruitment precedes synchronization. The model generates testable predictions for future clinical and experimental studies, provides a numerical dissection of how mesoscopic LFP features map onto microscopic neuronal dynamics, and serves as a computational building block for future circuit-level models that can guide brain stimulation strategies tailored to the patient-specific dynamical state of the STN. Author summaryIn Parkinsons disease, local field potentials (LFP) from the subthalamic nucleus (STN) contain two prominent rhythms: a slow beta rhythm (13-30 Hz) and fast oscillations (200-400 Hz). In the parkinsonian state, these rhythms become coupled, with fast oscillation amplitude varying systematically with beta phase, a relationship absent in the medicated state. We built a biophysical spiking neuron network model that captures two key effects of dopamine depletion on STN neuronal activity: changes in the intrinsic neuronal excitability and changes in synaptic coupling strength. The model produces fast oscillations from rapid intraburst firing, while the slow beta rhythm and its coupling to fast oscillations emerge with the onset of synchronized bursting across the population. Importantly, the frequency of the fast oscillations shifts continuously depending on both parameters, explaining a puzzling clinical observation that these oscillations change frequency between medication states. The model also reproduces the modulation pattern in the spike-triggered average of HFO envelope amplitude reported in patient recordings, confirming consistency with single-unit observations as well as LFP-level spectral features. By mapping how multi-timescale LFP spectral features relate to the dynamical regime of the underlying neuronal population, this work offers a framework for brain stimulation strategies informed by patient-specific dynamical states.